Novel symbionts and potential human pathogens excavated from argasid tick microbiomes that are shaped by dual or single symbiosis

Graphical abstract


Introduction
Ticks are obligatory blood-feeding arachnids and the second major vector transmitting pathogens to humans and animals [1]. Studies have shown that ticks pose a threat to human and animals not only through transmitting a wide range of pathogens but also https by causing blood loss, allergic reactions, and paralysis [2][3][4]. There are three tick families: Ixodidae (ixodid or hard-bodied ticks), Argasidae (argasid or soft-bodied ticks), and Nuttalliellidae (only one species) [5]. While each stage of female ixodid ticks feed only once and may have a diverse habitat, male ixodid ticks may feed several times. Some stages of argasids feed for a few hours or not at all (for example, Ornithodoros moubata larvae molt into nymphs without blood meals) [6] and live in hidden sites near their hosts [7]. More than 900 tick species have been reported [8] with more than 186 species classified as argasid species [9].
Generally, ticks harbor a wide variety of microbes, including endosymbionts, commensals, and tick-borne pathogens (TBPs) [10] that represent a complex microbiome. Some endosymbionts such as Coxiella, Francisella, and Rickettsia species are essential for ticks to provide the nutrients lacking in their blood meals [11]. In most cases, these symbionts are maternally inherited and have been associated with their vector hosts for millions of years [12]. Recently, studies describe more than seven endosymbiont bacterial genera from several tick species [13]. In addition, ticks harbor a wide variety of commensal bacteria from the surrounding environment and the host skin [14] such as Bacillus, Mycobacterium, Methylobacterium, and Pseudomonas [15]. Ticks are also frequently co-infected with several TBPs that can cause diseases in humans and animals [16]. Although many studies have described the microbiome of several tick species [17][18][19], only few have discussed the complexity of the epidemiology and network of interactions of these microbes within ticks [20][21][22][23][24]. For example, significant ecological associations have been detected between the nonpathogenic microbes and pathogens in some ixodid ticks [25,26]. To solve this complex network, an integration of metataxonomic approaches, bioinformatics and experimental design is required. Hence, an investigation on the rarely described microbiome and microbial interactions in argasid ticks is timely required to understand the role they play in the transmission of pathogens.
The global number of newly described argasid ticks has been increasing [27][28][29]. These argasid ticks have a high opportunity to acquire and transmit pathogens to their hosts in each developmental stage [14]. For instance, argasid ticks harbor several Rickettsia spp. [17,[30][31][32], Borrelia spp. [33][34][35], and Ehrlichia spp. [36]. To date, most of the large-scale studies on tick microbiomes have focused only on ixodid ticks [37] and little has been discussed about argasid ticks. Unlike ixodid ticks, ticks of the family Argasidae are nidicolous, having long life spans (6-9 years), and more frequent and fast feeders (several minutes to hours) [38][39][40]. Sequentially, all life stages of these ticks could coexist in one habitat, which can lead to maintaining several bacterial species throughout the developmental stages and sexes [41].
In Japan, a total of four argasid ticks have been described as Argas japonicus, Carios vespertilionis, Ornithodoros capensis, and Ornithodoros sawaii [42]. However, little is known about the role of these ticks in the transmission of TBPs to both humans and animals in a country where many human cases of tick-borne diseases (TBDs) are reported [43]. In addition, characterization of the diversity and composition of microbes associated with several argasid tick species is required to understand the evolutionary, ecological, and physiological functions of the tick-associated microbiome [44]. In addition, the study of tick microbial communities can help to improve the control strategies of ticks and TBPs [45] and identify novel bacterial species that may be pathogenic to humans and animals [46]. Recently, a study described the microbial diversity in A. japonicus samples collected from China [17], but the dominant endosymbiont bacteria was not clearly identified in these samples. In Mexico, two bacteria (Midichloria and Coxiella) were suggested to be the endosymbionts of O. turicata [47] based on the analysis of 17 adult ticks collected from the Bolson tortoise (Gopherus flavo-marginatus). Similarly, the microbiomes of seabird ticks (O. maritimus and O. capensis) were found to be dominated by Coxiella endosymbionts [41,48]. However, the microbiome of O. moubata, which is considered the most studied among argasid ticks, is dominated by two endosymbionts, Coxiella and Francisella [49,50].
In this study, four argasid tick species, A. japonicus, C. vespertilionis, O. capensis, and O. sawaii were used to characterize the microbiome composition and diversity in argasid ticks from Japan. To understand the relationship between the tick microbiome and the surrounding environment, we investigated the microbiome of the African O. moubata from a laboratory colony maintained at the University of Tsukuba, Japan. We report isolation and genetic characterization of a potentially novel endosymbiont detected in A. japonicus. Further, we targeted flagellin protein gene (flaB), citrate synthase gene (gltA), and chaperonin gene (groEL) to characterize a Borrelia and a potential novel Ehrlichia species that we detected in our samples. Collectively, our study can help to better understand the interrelationships between argasid ticks, their associated microbiomes, and the identification of potentially pathogenic bacteria.
All ticks were processed alive under sterile conditions. Tick DNA was extracted as described elsewhere [53]. In brief, the tick surface was washed once with 70% ethanol and twice with sterile PBS.

Bioinformatics processing
We used the Quantitative Insights Into Microbial Ecology 2 Software (QIIME2) (version 2020.2) [56] as a platform for microbiome bioinformatics in argasid ticks. The raw sequencing data were obtained from BaseSpace (Illumina), and then demultiplexed, quality-checked, and filtered using q2-demux plugin followed by denoising with DADA2 pipeline (version: 2019.10) [57]. The obtained amplicon sequence variants (ASVs) were aligned using q2-alignment plugin by mafft [58]. A phylogeny was constructed using q2-phylogeny plugin by fasttree2 [59]. We selected a sampling depth of 7,760 reads for comparing the diversity analysis using q2-diversity plugin in QIIME2 between the examined argasid tick species and excluded one sample from O. sawaii due to low number of reads. In addition, we compared the diversity values according to sex and stage within each argasid tick species using a sampling depth of 27342, 9442, 17,878, 14,069, and 7,760 for A. japonicus, C. vespertilionis, O. capensis, O. sawaii, and O. moubata, respectively. We calculated alpha diversity based on Shannon diversity [60] (accounts for richness and evenness), Faith's Phylogenetic Diversity (Faith's PD) [61] (accounts for the branch lengths of the rooted phylogenetic tree connecting the detected bacterial species), observed features [62] (accounts for the number of taxonomic groups), and Pielou's evenness [63] (accounts for the abundance of species relative to other species in a given community). We exported and visualized the results in R using the qiime2R, ggplot2, and phyloseq packages [64]. In addition, beta diversity was calculated based on unweighted UniFrac distance [65] (presence or absence of observed species), weighted UniFrac distance [66] (abundance of observed species), Jaccard similarity index [67] (the number of shared species divided by the total number of species detected in the sample sets), and Bray-Curtis dissimilarity [68] (differences between sites in terms of species counts present in those sites) analyses using QIIME2. Furthermore, visualization of clustering of ASVs according to species and sex/stage was performed by a Principal Coordinates Analysis (PCoA) using EMPeror plugin in QIIME2 [69] and R as described above. The taxonomy was assigned using q2-feature-classifier plugin [70] with classify-sklearn naïve Bayes taxonomy classifier and SILVA classifier reference sequences (release 132). The frequency method in the Decontam package [71] in R (version 4.0.2) [72] was implemented with a threshold of 0.1 to identify the likely contaminants introduced during processing. Finally, ASVs identified as archaea, eukaryota, and potential contaminants as well as sequences not assigned to kingdom level were removed manually for further analysis in QIIME2 using sequence identifiers. A heatmap phylogenetic tree was constructed using the heatmap method in QIIME2 [73]. We visualized the differential abundance of the 30 most abundant taxonomic groups using taxa_heatmap function in the qiime2R package in R (version 2.13.0). Finally, we identified the bacteria contributing to the dissimilarity of the microbiome among the tick groups by implementing the linear discriminant analysis effect size (LEfSe) in the Huttenhower lab Galaxy pipeline [74]. LEfSe was implemented to compare among (between different species for the same sex and stage) and within (between different sex and stage for each species) tick species variations.

Statistical analyses
We estimated the statistical differences in alpha diversities among argasid tick species using a generalized linear model (GLM). The response variable was alpha diversity including Shannon diversity calculated as the exponential function of Shannon entropy [75], Faith's PD, observed features, and Pielou's evenness with tick species and sex/stage (female, male, and nymph) as fixed effect variables. We checked various model assumptions, including normality of residuals, normality of random effects, linear relationship, homogeneity of variance, and multicollinearity, using the check_model function in the performance package (version 0.8.0) in R [76]. We used the emmeans package in R [77] as a post hoc method to show the resulted pairwise comparisons. The modeling was performed using the stats package version 4.0.3 in R [72]. We then analyzed the effect of species and sex/stage on beta diversity among all argasid ticks using Adonis Permutational multivariate analysis of variance (PERMANOVA) with 999 permutations. Afterwards, we ran a series of pairwise comparisons between all argasid species. We also tested the significance in beta diversity differences within each argasid tick species according to sex/stage using pairwise PERMANOVA [78] with 999 permutations.

Isolation and characterization of potential endosymbionts
A total of 30 tick homogenates (5 C. vespertilionis, 12 A. japonicus, 6 O. capensis, 5 O. sawaii, and 2 O. moubata) were selected to isolate potential endosymbionts using arthropod cell lines as described previously [53]. In brief, tick homogenates (5 lL each) were inoculated to ISE6 (Ixodes scapularis embryos) cell and C6/36 (Aedes albopictus larvae) cell cultures seeded in 24-well culture plates. The culture medium was changed every 3 and 7 days for C6/36 and ISE6 cells, respectively. Every 2 weeks, one hundred lL of culture suspension was blindly passaged into new wells containing uninfected cells. The experiment was terminated at 8 weeks post inoculation. When fungal or bacterial contamination was observed, the contaminated wells were sterilized with 10% hypochlorous acid.
In addition to microscopic observations to detect cytopathic effects caused by bacterial infections, crude DNA extracted from 10 lL of culture suspension at 4-and 8-weeks post inoculation using heat alkali treatment was tested for bacterial infection by universal PCR amplifying bacterial 16S rRNA gene. PCR was carried out using the high-fidelity KOD-Plus-Neo DNA polymerase (Toyobo, Osaka, Japan) in a 25.0 lL reaction mixture containing 2.5 lL of 10x Buffer for KOD-Plus-Neo, 300 nM of each fD1 and Rp2 primers [79], 2.5 lL of dNTPs, 1.5 lL of 25 mM MgSO 4 , 0.5 unit of KOD-Plus-Neo DNA polymerase, 2.0 lL of template DNA and molecular grade water was used to adjust the volume. The reaction conditions were set at 94°C for 2 min and 30 cycles of 98°C for 10 sec, 55°C for 30 sec and 68°C for 30 sec, and final extension at 68°C for 2 min (Table 1).
Borrelia-positive argasid tick samples were examined by nested PCR targeting a 345 bp segment of Borrelia spp. flaB using the primer pairs BflaPAD and BflaPDU for the first PCR and BflaPBU and BflaPCR (Table 1) for the secondary PCR [81]. PCR cycling conditions were implemented as previously described [82].
Finally, we examined one Wolbachia-positive argasid tick sample by PCR targeting a 1,395 bp segment of Wolbachia spp. 16S rRNA gene using the primer pair EHR16SD [83] and 1513R [79] ( Table 1). The PCR products were electrophoresed in a 1.5% agarose gel stained with Gel-Red and visualized under UV light.

Sanger sequencing and phylogenetic analysis
The resulting PCR products were purified using the NucleoSpin Gel and PCR Clean Up Kit (Takara Bio Inc.) or ExoSAP-IT PCR Product Cleanup Reagent (Thermo Fisher Scientific, Waltham, MA, USA) and sequenced by using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) and an ABI Prism 3130 Â genetic analyser (Applied Biosystems) according to the manufacturers' instructions.
The obtained sequences were imported to Geneious v10.2.6 (Biomatters Ltd., Auckland, New Zealand) and the primer sequences were removed. The assembled sequences were submitted to the DNA Data Bank of Japan (DDBJ) under the accession numbers LC649932-LC649949.
The nucleotide Basic Local Alignment Search Tool (BLASTn) (accessed on 5 June 2021, https://blast.ncbi.nlm.nih.gov/Blast.cgi) was used to compare the sequences in this study with public databases. We used MAFFT v7.450 software [58] to align the sequences with the previously published related ones. The best fit model for the phylogenetic analysis was estimated using MEGA X software [84] and the trees were built using PHYML v3.3 software [85] with the branches supported by Bootstrap tests.

Laboratory colonies of Ornithodoros moubata harbor less diverse microbiome than field collected argasid ticks
A total of 8,394,800 raw paired-end reads were generated by the Illumina MiSeq from a total of 137 argasid tick and eight negative control samples. The reads were truncated and trimmed according to the demultiplexing summary. The DADA2 quality control analysis and Decontam package output resulted in 4,536,908 high quality paired-end reads classified into 4,279 features. The maximum and minimum frequencies per sample were 61,825 and 4,260, respectively. The lowest total number of reads and features was detected from O. moubata (Table 2, Fig. S1).

Tick species and sex/stage elucidate variations in the associated microbiome diversity
There were no potential issues in our GLMs based on the visualized values of model assumptions that were used to test significance in alpha diversity among the tested argasid ticks (Fig. S2-5). The alpha diversity of tick microbiome among the examined argasid ticks using Shannon, Faith's PD, observed features, and Pielou's evenness analyses is significantly affected by the tick species variation as tested by GLM (p < 0.01) ( Fig. 1 and Table S2-5). Interestingly, we found that the four alpha diversity metrices were significantly lower in Argas and Carios species when compared to Ornithodoros spp. collected from the environments (p < 0.01) (Fig. 1, Fig. S6, and Table S2-5) except that Shannon diversity results of A. japonicus and C. vespertilionis were not significantly different from O. moubata and O. sawii, respectively. In addition, Pielou's evenness of C. vespertilionis was not significantly different from O. moubata. We also found that all alpha diversity metrices, except the observed features, of O. capensis were not significantly different from O. sawaii. The analysis also showed significant differences in Shannon diversity, Faith's PD, and the observed features of adult argasid ticks when compared to nymphs.
Four alpha diversity metrices showed non-significant correlation with the sex and stage variations within all examined argasid species except for O. capensis where the values in nymphs were significantly lower than in adult stages (Table S6, Fig. 1, and Fig. S6). The estimated Pielou's evenness diversity within O. moubata, the only laboratory colony in this study, were significantly lower in the nymphal stage than in the adult females and males (p 0.01) (Table S6 and Fig. S6). In addition, adonis PERMANOVA showed that the calculated beta diversity metrices including the unweighted UniFrac, weighted UniFrac, Jaccard, and Bray-Curtis analyses were significantly variable according to argasid species and sex/stage (p < 0.01) (Table S7). Pairwise PERMANOVA showed that the four beta diversity metrices were significantly different between all argasid species (Table S8 and Fig. S7) and the effect of sex/stage on beta diversity within each species was not consistent as explained among different species. Specifically, only Jaccard was significantly affected by sex and stage variation in A. japonicus (p < 0.01) (Table S9) and C. vespertilionis (p < 0.05) (Table S10), respectively. The pairwise analysis of O. capensis showed signifi-  (Table S13). Our PCoA plots confirmed the pairwise PERMANOVA results as samples from each argasid species clustered together separate from the others (Fig. 1). In addition, there was an overlap in the clustering of sex and stage within species (Fig. 1 and Fig. S8).

Dual symbiosis in argasid ticks and the detection of novel endosymbionts in Argas japonicus
The most abundant four bacterial families were completely different among the examined argasid tick species (Table 3 and Fig. S9). We presented the total number of reads for each bacterial feature detected in this study in Table S14. According to the list of the most common maternally inherited bacteria in ticks [44], the microbiome of each examined argasid species included at least one potential endosymbiont and three species (A. japonicus, C. vespertilionis, and O. moubata) had two potential endosymbionts ( Fig. 2 and Fig. S10-14).
The relative abundance and LEfSe analyses of the bacterial genera showed a unique presentation of the potential endosymbionts, commensals, and pathogens in the examined five argasid species. A bacterial genus belonging to the family Piscirickettsiaceae (identified as Candidatus Endoecteinascidia in SILVA classifier) and Occidentia species were the most abundant in A. japonicus (mean relative abundance = 53.0% and 27.2%, respectively) ( Table 4). However, the taxa bar plot showed a negative correlation between both species in A. japonicus (Fig. 2 and Fig. S10). In addition, A. japonicus samples were divided into two clusters according to the relative abundance of the potential pathogenic or endosymbiotic bacterial taxa, where each cluster was dominated by either  Candidatus Endoecteinascidia or Occidentia species (Fig. 3). The most abundant bacterial genera in C. vespertilionis were Rickettsia and Coxiella (mean relative abundance = 64.0% and 21.4%, respectively) ( Table 4, Fig. 3, and Fig. S11).
Although the microbiome of O. capensis was mainly dominated by commensal bacteria (Fig. S12), we found that the mean relative abundance of the potential endosymbiont Coxiella sp. was 2.6% and detected in 93.9% (46/49) of the examined O. capensis samples (Table 4). Similarly, commensals dominated the microbiome of O. sawaii samples (Fig. S13) and the mean relative abundance of Coxiella sp. was 7.2% detected from 100% (20/20) of the examined O. sawaii samples (Table 4). In contrast, the microbiome of the laboratory colony of O. moubata was dominated by the bacterial genera Francisella and Coxiella (Table 4 and Fig. S14) (mean relative abundance = 64.2% and 34.9%, respectively).
LEfSe analysis detected differences in the community composition between the examined argasid species for the same sex and stage ( Fig. S15-17). Far fewer discriminant bacterial taxa were detected within each argasid species when compared according to the sex and life stage. Although the majority of these bacterial taxa were commensals, we detected significant differential abundances for endosymbionts and pathogens in relation to sex and stage ( Fig. 3 and Fig. S18-19). Notably, the relative abundance of Occidentia sp. was higher in female than male A. japonicus. Female O. capensis showed higher abundance of Coxiella sp. than in males from the same species. Interestingly, the abundance of Francisella sp. in the nymphal stage was significantly higher than in both male and female O. moubata. In contrast, adult O. moubata had higher abundance of Coxiella than nymphs.

Isolation and molecular characterization of Occidentia massiliensis from Argas japonicus
During an 8-week incubation period, a bacterial isolate was obtained from one single A. japonicus sample inoculated into ISE6 cells (Fig. 4). A total of 14 and 27 samples were contaminated with environmental bacteria/fungi for ISE6 and C6/36 cells, respectively. The sequencing analysis of 16S rRNA gene PCR amplicon of the isolate (GenBank accession number: LC649932) showed 99% identity to Oc. massiliensis strain Os18 (GenBank accession number: NR_149220), which was isolated from Ornithodoros sonari from Senegal. The phylogenetic tree showed that the obtained sequence clustered with Oc. massiliensis and formed a separate clade between the previously published Orientia and Rickettsia species (Fig. 5a).

Molecular characterization of a potential novel endosymbiont belonging to the order Thiotrichales
Although this bacterium was identified as Candidatus Endoecteinascidia by the SILVA classifier in our taxonomy analysis, the obtained nearly full 16S rRNA gene sequences were identical to   each other and showed 96% identity to uncultured Thiotrichales (GenBank accession number: LC278460) isolated from Asterias amurensis from Hokkaido, Japan. The phylogenetic tree showed that the obtained four sequences (GenBank accession numbers: LC649933-LC649936) clustered with the previously published sequences from uncultured bacteria belonging to the order Thiotrichales and formed a clade adjacent to the cluster of Francisella species (Fig. 5b).

Molecular characterization of potential novel pathogens in Ornithodoros species
The partial gltA gene sequences (GenBank accession numbers: LC649937-LC649941) of the detected Ehrlichia sp. (n = 5) were identical and showed 85.2% identity to E. canis (GenBank accession number: AY647155) reported from dogs in Italy. The partial groEL gene sequences (GenBank accession numbers: LC649942-  The phylogenetic analysis showed that the detected Borrelia sp. sequences clustered together with Borrelia turicatae group (Fig. 5e). The obtained nearly full 16S rRNA gene sequence of Wolbachia sp. (GenBank accession number: LC677220) showed 99.4% identity to Wolbachia endosymbiont of Folsomia candida (GenBank accession number: EU330430) isolated from USA. The phylogenetic tree showed that the obtained sequence clustered with the previously published sequences from a Wolbachia endosymbiont of other arthropods (Fig. 5f).

Discussion
During the past decade, the research on vector-associated microbiomes has been largely expanding due to improvements of sequencing technology and bioinformatic tools in this field that has contributed greatly to better understanding of tick systematics [86,87] and the emergence of novel TBPs [88]. However, little is known about microbiomes of argasid ticks and the drivers affecting its composition [41]. In our study, microbiomes of 137 argasid ticks of 5 different species were investigated and the association with species and sex/stage variations in Japan was explained. Two novel putative endosymbionts, together with novel Ehrlichia, Borrelia, and Wolbachia were molecularly characterized among these argasid samples.
Although our study did not include field collected O. moubata samples, the examined laboratory colony showed significantly lower number of features as compared to those in the field collected samples from the other four species. The difference in microbiome between field and laboratory collected ticks has previously been demonstrated in Amblyomma maculatum [89], where the microbiome diversity was higher in the field collected A. maculatum than in the laboratory reared ones. In another study, laboratory colonization of Dermacentor andersoni resulted in stabilization of the associated microbiome [90]. Hence, our results support the hypothesis that laboratory colonization of ticks can reduce the role of ecological factors in changing the associated microbiome in a way that can improve the biological control of ticks by targeting the essential remaining endosymbionts after colonization. This strategy has helped in the production of germ-free mosquitoes and the acceleration of the study of microbiome in insects [91].
Few studies have demonstrated the effect of species variation on argasid tick microbiome [10,22,92] due to diverse animal hosts and geographical distribution. We have showed that the differences in microbiome composition and richness in this study were best explained by variations in tick species. This strong speciesspecific variation in the microbiome of argasid ticks can be attributed to differences in collection sites [93,94] as well as potential animal hosts [95]. In this study, samples representing each argasid species were collected from a specific collection site, which hindered testing the effect of spatial variation within the same species. Recorded host preferences are variable among the argasid species employed in this study. Briefly, the host animals for A. japonicus, C. vespertilionis, O. capensis, and O. sawaii are the common house martin (Delichon urbicum), the Japanese short-tailed bat (Ep. japonensis), the wedge-tailed shearwater (Puffinus pacificus), and the Swinhoe's petrel (Oce. monorhis) [96,97], respectively. This animal host variation may explain the significant differences in microbiome composition and richness between the examined species of ticks. In general, argasid ticks feed several times from their animal hosts during their life cycles [98]. Hence, many commensal bacteria and pathogens can be transmitted to these ticks from the skin and blood of their animal hosts during feeding [13,14]. The taxonomy composition analysis showed that Coxiella was the potential endosymbiont for 4 out of 5 argasid species. Only A. japonicus was dominated either by Oc. massiliensis (11/38) or an uncultured Thiotrichales bacterium (27/38). In addition, dual symbiosis was observed between Rickettsia and Coxiella or Francisella and Coxiella in C. vespertilionis and O. moubata, respectively. Some of the members of Coxiella, Rickettsia, and Francisella (obligate intracellular bacteria transmitted vertically) have a major nutritional role in ticks and have been detected from ixodid and argasid ticks [13,47,[99][100][101][102]. Previously, high dual abundance of Rickettsia and Coxiella was detected in O. capensis and Ornithodoros maritimus [41,48]. Furthermore, experimental elimination of Francisella endosymbiont in O. moubata has prevented the development of adult females [101]. Differences in microbiome were not only limited to between different argasid species but was also detected among individual ticks belonging to the same species (Fig. 2 and Fig. S10-14). These differences may be due to extrinsic or intrinsic variables such as the tick general physiological health, blood meal sources and histories, presence or absence of pathogenic bacteria and protozoa, tick life stage and sex, and tick genetic backgrounds [26,[103][104][105]. Our study design cannot illustrate all these factors to explain this variability among individuals from the same species, but we showed that infection with potential pathogens, sex, and stage can contribute to this pattern of microbiome variability.
Generally, co-occurrence of endosymbionts in one host may have several potential beneficiary consequences such as nutritional and protective functions [106]. For example, Midichloria could compensate the biotin biosynthesis pathways that are deficient in the co-symbiont Francisella in Hyalomma marginatum [107]. Another example is the agonist interaction between two endosymbionts through synthesizing different important nutrients in the glassy-winged sharpshooter (Homalodisca coagulata) [108]. In Aphidius ervi, the primary symbiont, Buchnera aphidicola, contributes mainly to the nutritional needs of the host, whereas the secondary symbiont, 'Candidatus Hamiltonella defensa', encodes toxins that protect the aphid from parasitoids [109]. In ixodid ticks, Rickettsia, Coxiella, and Francisella endosymbionts were found to co-exist in Dermacentor reticulatus and Ixodes ricinus collected in Slovakia [110]. Furthermore, a comparison between Coxiella and Francisella endosymbionts from ixodid ticks has revealed that their metabolism is highly similar, but only Francisella can produce more products than Coxiella including cysteine, threonine, tyrosine, tryptophan, phenylalanine, serine, asparagine, glutamine, proline, and heme, while Coxiella generates thiamine [111]. We assume that bacterial endosymbionts in argasid ticks are highly related to the blood meal source, especially that the metabolites provided are highly variable between the nucleated red blood cells (RBCs) in birds and reptiles, and the non-nucleated RBCs in mammals [112]. This can be supported by our results that C. vespertilionis and O. moubata ticks, which feed on blood meals containing mostly non-nucleated cells from Japanese short-tailed bats and mammals, were associated with dual symbiosis while the remaining three argasid ticks, which feed on blood meals containing mostly nucleated cells from birds, were associated with a single potential bacterial symbiont. Future investigations are required to confirm the effect of the blood meal source on the evolution of bacterial symbiosis in argasid ticks.
Previously, Oc. massiliensis was isolated from O. sonrai collected from rodent burrows in Senegal [113] and from Africaniella transversale collected from Python regius imported from Senegal to the United Arab Emirates [114], but no clear information was provided on its role as an endosymbiont. Because O. sonrai and Af. transversale are known to infest reptilian hosts during their life cycles [115,116], it appears that pythons and rodents (pythons' main prey) are equally susceptible to Oc. massiliensis infections. However, our study has isolated and characterized Oc. massiliensis for the first time from A. japonicus samples collected from Japan, suggesting a wider geographical distribution and host range of this endosymbiont than previously expected. We showed that 16S rRNA gene sequences from the isolated Oc. massiliensis is closely related to Rickettsia endosymbionts of the silverleaf whitefly (Bemisia tabaci) [117] and the Indian grain aphid (Sitobion miscanthi) [118,119], which supports the possibility that this bacterium can have a role as an endosymbiont in A. japonicus. Although we were not successful in isolating the other potential endosymbionts of A. japonicus that belong to the order Thiotrichales, we showed that this uncultured bacterium 16S rRNA gene sequence is closely related to the clade containing Francisella spp. that are considered as endosymbionts of several Argas species [120].
A previous study showed that the microbiome of female Ixodes ticks was less diverse than in males [94] due to the high abundance of Rickettsia species in females [80]. Although we have not included larvae in our investigation, we detected significant correlations between the life stage and the microbiome diversity in O. capensis and O. moubata. The explanation of this correlation may depend on the balance between the acquired environmental bacteria and the burden of endosymbionts within argasid ticks. This can explain the inconsistent alpha diversity values between different life stages in the tick species of this study ( Fig. 1 and S6). We showed that the nymphal stage of O. capensis and O. moubata presented lower alpha diversity values than adults, but no significant correlations were detected within the remaining argasids according to sex and stage variations. There are contrasting opinions regarding the interpretation of the correlation between life stage and microbiome diversity in ticks. One opinion is that there are positive correlations between alpha diversity and the tick life stage [121], except that females have lower diversity due to the colonization by endosymbionts such as Rickettsia species that lead to competition with other bacteria for resources [22,122]. Another opinion is that alpha diversity decreases by tick life stage due to the dilution effect through bloodmeal feeding in I. pacificus [123]. Our results from four field collected and one laboratory colony of argasid ticks suggest that it is difficult to find a consistent trend of microbiome diversity in relation to tick life stage that fits all tick species, and the final interpretation should be related to the life history of the samples and the balance between environmental bacterial exposure and endosymbiont ecology.
The LEfSe analysis showed that adult females of O. capensis and O. moubata have significantly higher abundance of Coxiella. Interestingly, the relative abundance of Francisella in O. moubata was significantly higher in nymphs than in adults. The relationship between Francisella and Coxiella in O. moubata with the life stage is not fully covered but it was suggested that Francisella has replaced Coxiella in some ticks [44,124]. A similar trend was observed in Amblyomma americanum where males harbored higher abundance of Rickettsia than females wihch had higher levels of Coxiella endosymbionts [125]. This may indicate that males have higher numbers of Rickettsia to facilitate locomotion to increase the mating opportunities [126]. However, it is not clear yet if Francisella has any role in locomotion facilitation in argasids. Although we detected significant variations in microbial diversity and richness among several developmental stages and sexes within the same argasid tick species, our analysis revealed that some bacterial dynamics are solid throughout the tick life stage. Indeed, argasid ticks are nidicolous but can disperse among several colonies of seabird species during the movement of their hosts [127], where new bacterial species can be acquired from the different environment [41].
Most of the previously published studies on Ehrlichia species have focused on ixodid ticks and animals [128,129]. We detected Ehrlichia from O. capensis and O. sawaii which showed 85 -91% identity to the previously published gltA and groEL gene sequences of Ehrlichia species. The number of the reports on uncultured Ehrlichia spp. has been increasing worldwide. Recently, a study has characterized a total of 11 newly identified Ehrlichia genotypes detected in Haemaphysalis ticks in Japan [130]. To the best of our knowledge, this is the first study to report Ehrlichia spp. from argasid ticks in Japan. A closely related species to Candidatus Ehrlichia khabarensis has been detected from C. vespertilionis in the UK [36] and a new species related to E. canis group (Ehrlichia sp. AvBat) was detected from the same tick species from France [131]. Due to the relative lack of available data regarding the prevalence of Ehrlichia spp. in Ornithodoros ticks, more investigations are needed in the future to identify the role of argasid ticks as vectors for TBPs and the pathogenicity of the detected Ehrlichia sp. to humans and animals.
Several bacterial TBDs, such as Japanese spotted fever, Lyme disease, relapsing fever, and scrub typhus, are endemic in Japan [132]. We detected a Borrelia sp. in O. capensis that clustered with the relapsing fever Borrelia group including Borrelia parkeri, B. turicatae, Borrelia venezuelensis, and Borrelia johnsonii. Previous studies in Japan have detected several Lyme disease borreliae (eg, Borrelia japonica) and relapsing fever borreliae (eg, Borrelia miyamotoi) in ixodid ticks [82,133,134]. In addition, Borrelia sp. K64 has an identical partial flaB gene sequence with the one obtained in our study, and it was reported from O. sawaii in Japan [135]. This Borrelia sp. has been detected recently from O. sawaii in South Korea [97], but our study is the first to detect this species in O. capensis. These findings support the possibility that migratory seabirds are the vertebrate hosts for this Borrelia species. Additionally, the presence of Borrelia sp. closely related to B. turicatae was highlighted in O. capensis from Algeria, which suggests that this Borrelia sp. is main-tained by seabirds and transmitted by the associated argasid ticks [136]. Hence, more investigations are required to know if this Borrelia sp. could represent a public health threat in Japan and South Korea.
Wolbachia species are highly diverse in invertebrates and represent an important group of endosymbionts in nematodes and insects [137][138][139]. Our study is the first to detect a Wolbachia sp. in O. sawaii. Previously, several Wolbachia spp. were detected from I. ricinus from France, Germany, Italy, and The Netherlands [140] and it was suggested that this bacterium could have negative impact on the tick microbiome due to the induction of immune reactions [141]. The presence of Wolbachia sp. in one of our O. sawaii samples can be due to the association of this argasid sample with an endoparasitoid as was shown for I. ricinus samples infected with Wolbachia sp. [142] due to the infestation by Ixodiphagus hookeri (Hymenoptera, Chalcidoidea, Encyrtidae). Although we did not investigate the presence of endoparasitoids in our samples, the phylogenetic analysis showed that the detected Wolbachia sp. in our study is closely related to the endosymbiont of I. hookeri, which supports the hypothesis that this argasid tick sample was probably infested with an endoparasitoid species that uses this Wolbachia sp. as an endosymbiont.

Conclusions
By targeting the bacterial communities in 137 argasid ticks collected from multiple locations in Japan, we were able to provide insights on the microbiome of four natural and one laboratory reared argasid species. We showed that the microbial community structures within each argasid species is distinct from the other and can be affected by life stage variation. Furthermore, our study illustrated that only a few bacterial genera differ in abundance within the same argasid species for different life stages. Specifically, Occidentia, Coxiella, and Francisella endosymbionts showed significant correlations with the argasid tick life stage variations. We documented the isolation and molecular characterization of novel endosymbionts in A. japonicus, including Oc. massiliensis and a bacterium belonging to the order Thiotrichales. The two newly identified potential symbionts expand the current identified phylogenetic clades of tick endosymbionts. We detected and characterized a potentially novel Ehrlichia sp. from O. capensis and O. sawaii. Additionally, we detected a relapsing fever Borrelia sp. for the first time from O. capensis. Our study should encourage future investigations on argasid ticks to enrich our understanding about microbiome ecology and possible manipulations to provide next generation control strategies for ticks. Likewise, investigations are needed to identify the reservoir capacity of argasid ticks to TBPs with zoonotic and veterinary importance, considering proper isolation and characterization of the detected pathogens.

Availability of data and materials
The raw sequence data were deposited in the DNA Data Bank of the Japan Sequence Read Archive under the DRA accession number: DRA013238.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.